#4 Генетика ~ Продолжительность жизни
rm(list=ls(all=TRUE))
#wd = getwd()
#wd = paste(wd, '/mtdna-mammalian-evolution/Body/2Derived',sep='')
#setwd(wd)
library(ggplot2)
data = read.table("../../Body/2Derived/CytB_with_ecology.csv", sep='\t', header = TRUE)
###### КОРЕЛЯЦИИ
#(i) FractionOfSyn ~ -GenerLength;
cor.test(x = data$GenerationLength_d, y = data$FractionOfSyn , method = "spearm")
## Warning in cor.test.default(x = data$GenerationLength_d, y =
## data$FractionOfSyn, : Cannot compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: data$GenerationLength_d and data$FractionOfSyn
## S = 10107033, p-value = 1.771e-13
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.3786509
GenerationLength_FractionOfSyn_fit <- cor.test(x = data$GenerationLength_d, y = data$FractionOfSyn)
ggplot(data, aes(x = log(GenerationLength_d), y = FractionOfSyn ,col = factor(Order) ))+
geom_point(size = 2)

ggplot(data, aes(x = log(GenerationLength_d), y = FractionOfSyn , fill = Order)) +
geom_boxplot()

#(ii) FractionOfNonsyn ~ +GenerLength;
cor.test(x = data$GenerationLength_d, y = data$FractionOfNonsyn , method = "spearm")
## Warning in cor.test.default(x = data$GenerationLength_d, y =
## data$FractionOfNonsyn, : Cannot compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: data$GenerationLength_d and data$FractionOfNonsyn
## S = 4555175, p-value = 1.771e-13
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.3786509
GenerationLength_FractionOfNonsyn_fit <- cor.test(x = data$GenerationLength_d, y = data$FractionOfNonsyn)
ggplot(data, aes(x = log(GenerationLength_d), y = FractionOfNonsyn ,col = factor(Order) ))+
geom_point(size = 2)

ggplot(data, aes(x = log(GenerationLength_d), y = FractionOfNonsyn , fill = Order)) +
geom_boxplot()

#(iii) SummOfAllGrantham ~ +GenerLength;
cor.test(x = data$GenerationLength_d, y = data$SummOfAllGrantham , method = "spearm")
## Warning in cor.test.default(x = data$GenerationLength_d, y =
## data$SummOfAllGrantham, : Cannot compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: data$GenerationLength_d and data$SummOfAllGrantham
## S = 6351234, p-value = 0.01195
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.1336593
GenerationLength_SummOfAllGrantham_fit <- cor.test(x = data$GenerationLength_d, y = data$SummOfAllGrantham)
ggplot(data, aes(x = log(GenerationLength_d), y = SummOfAllGrantham ,col = factor(Order) ))+
geom_point(size = 2)

ggplot(data, aes(x = log(GenerationLength_d), y = SummOfAllGrantham , fill = Order)) +
geom_boxplot()

#(iv) MedianOfAllNonsyn ~ +GenerLength
cor.test(x = data$GenerationLength_d, y = data$MedianOfAllNonsyn , method = "spearm")
## Warning in cor.test.default(x = data$GenerationLength_d, y =
## data$MedianOfAllNonsyn, : Cannot compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: data$GenerationLength_d and data$MedianOfAllNonsyn
## S = 6992814, p-value = 0.3874
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.04614444
GenerationLength_MedianOfAllNonsyn_fit <- cor.test(x = data$GenerationLength_d, y = data$MedianOfAllNonsyn)
ggplot(data, aes(x = log(GenerationLength_d), y = MedianOfAllNonsyn ,col = factor(Order) ))+
geom_point(size = 2)

ggplot(data, aes(x = log(GenerationLength_d), y =MedianOfAllNonsyn , fill = Order)) +
geom_boxplot()

#(v) SummOfAllGrantham ~+SummOfAllNonsyn (проверка на вшивость)
cor.test(x = data$SummOfAllGrantham, y = data$MedianOfAllNonsyn, method = "spearm")
## Warning in cor.test.default(x = data$SummOfAllGrantham, y =
## data$MedianOfAllNonsyn, : Cannot compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: data$SummOfAllGrantham and data$MedianOfAllNonsyn
## S = 770600, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.8948862
GenerationLength_MedianOfAllNonsyn_fit <- cor.test(x = data$SummOfAllGrantham, y = data$MedianOfAllNonsyn)
ggplot(data, aes(x = log(SummOfAllGrantham), y = MedianOfAllNonsyn ,col = factor(Order) ))+
geom_point(size = 2)

ggplot(data, aes(x = log(SummOfAllGrantham), y = MedianOfAllNonsyn , fill = Order)) +
geom_boxplot()
## Warning: Removed 25 rows containing non-finite values (stat_boxplot).

#(vii) AverageGrantham ~ +GenerLength;
cor.test(x = data$GenerationLength_d, y = data$AverageGrantham , method = "spearm")
## Warning in cor.test.default(x = data$GenerationLength_d, y =
## data$AverageGrantham, : Cannot compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: data$GenerationLength_d and data$AverageGrantham
## S = 5006904, p-value = 1.106e-09
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.3170328
GenerationLength_AverageGrantham_fit <- cor.test(x = data$GenerationLength_d, y = data$AverageGrantham)
ggplot(data, aes(x = log(GenerationLength_d), y = AverageGrantham ,col = factor(Order) ))+
geom_point(size = 2)

ggplot(data, aes(x = log(GenerationLength_d), y = AverageGrantham, fill = Order)) +
geom_boxplot()

#(viii) KnKs ~ +GenerLength;
cor.test(x = data$GenerationLength_d, y = data$KnKs , method = "spearm")
## Warning in cor.test.default(x = data$GenerationLength_d, y = data$KnKs, : Cannot
## compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: data$GenerationLength_d and data$KnKs
## S = 4588296, p-value = 3.592e-13
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.374133
GenerationLength_KnKs_fit <- cor.test(x = data$GenerationLength_d, y = data$KnKs)
ggplot(data, aes(x = log(GenerationLength_d), y = KnKs ,col = factor(Order) ))+
geom_point(size = 2)

ggplot(data, aes(x = log(GenerationLength_d), y = KnKs, fill = Order)) +
geom_boxplot()

#(ix) DnDs ~ +GenerLength;
cor.test(x = data$GenerationLength_d, y = data$DnDs , method = "spearm")
## Warning in cor.test.default(x = data$GenerationLength_d, y = data$DnDs, : Cannot
## compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: data$GenerationLength_d and data$DnDs
## S = 4135683, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.4358718
GenerationLength_DnDs_fit <- cor.test(x = data$GenerationLength_d, y = data$DnDs)
ggplot(data, aes(x = log(GenerationLength_d), y = DnDs ,col = factor(Order) ))+
geom_point(size = 2)

ggplot(data, aes(x = log(GenerationLength_d), y = DnDs, fill = Order)) +
geom_boxplot()

##################### ПИКИ
library(ape)
library(geiger)
library(caper)
## Loading required package: MASS
## Loading required package: mvtnorm
library(stringr)
library(phytools)
## Loading required package: maps
library(picante)
## Loading required package: vegan
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-7
##
## Attaching package: 'vegan'
## The following object is masked from 'package:phytools':
##
## scores
## Loading required package: nlme
tree = read.tree("../../Body/1Raw/mammalia_species.nwk")
data$Species = sub(" ", "_", data$Species, ignore.case = FALSE, perl = FALSE, fixed = FALSE, useBytes = FALSE)
row.names(data) = data$Species
nrow(data) # 344
## [1] 353
tree_w = treedata(tree, data[, c('Species', 'AverageGrantham', 'KnKs','DnDs' ,'GenerationLength_d')],
sort=T, warnings=T)$phy
## Warning in treedata(tree, data[, c("Species", "AverageGrantham", "KnKs", : The following tips were not found in 'data' and were dropped from 'phy':
## Abeomelomys_sevia
## Abrocoma_bennettii
## Abrocoma_cinerea
## Abrothrix_andinus
## Abrothrix_hershkovitzi
## Abrothrix_illuteus
## Abrothrix_jelskii
## Abrothrix_lanosus
## Abrothrix_markhami
## Acerodon_celebensis
## Acerodon_jubatus
## Acinonyx_jubatus
## Acomys_chudeaui
## Acomys_cilicicus
## Acomys_cineraceus
## Acomys_ignitus
## Acomys_johannis
## Acomys_kempi
## Acomys_minous
## Acomys_nesiotes
## Acomys_ngurui
## Acomys_percivali
## Acomys_russatus
## Acomys_spinosissimus
## Acomys_subspinosus
## Acomys_wilsoni
## Aconaemys_fuscus
## Aconaemys_porteri
## Aconaemys_sagei
## Acrobates_pygmaeus
## Addax_nasomaculatus
## Aegialomys_xanthaeolus
## Aepeomys_lugens
## Aepyceros_melampus
## Aepyprymnus_rufescens
## Aeretes_melanopterus
## Aeromys_tephromelas
## Aethalops_aequalis
## Aethalops_alecto
## Aethomys_chrysophilus
## Aethomys_ineptus
## Aethomys_kaiseri
## Ailurops_ursinus
## Ailurus_fulgens
## Akodon_aerosus
## Akodon_affinis
## Akodon_albiventer
## Akodon_azarae
## Akodon_baliolus
## Akodon_boliviensis
## Akodon_budini
## Akodon_cursor
## Akodon_dayi
## Akodon_dolores
## Akodon_fumeus
## Akodon_iniscatus
## Akodon_juninensis
## Akodon_kofordi
## Akodon_latebricola
## Akodon_lindberghi
## Akodon_lutescens
## Akodon_mimus
## Akodon_molinae
## Akodon_mollis
## Akodon_mystax
## Akodon_orophilus
## Akodon_paranaensis
## Akodon_reigi
## Akodon_serrensis
## Akodon_siberiae
## Akodon_simulator
## Akodon_subfuscus
## Akodon_sylvanus
## Akodon_toba
## Akodon_torques
## Akodon_varius
## Alcelaphus_buselaphus
## Alcelaphus_caama
## Alcelaphus_lichtensteini
## Alces_alces
## Alces_americanus
## Alionycteris_paucidentata
## Allactaga_balikunica
## Allactaga_bullata
## Allactaga_elater
## Allactaga_euphratica
## Allactaga_firouzi
## Allactaga_hotsoni
## Allactaga_major
## Allactaga_sibirica
## Allactaga_sp._Kazakhstan
## Allactaga_toussi
## Allactaga_williamsi
## Allactodipus_bobrinskii
## Allenopithecus_nigroviridis
## Allocebus_trichotis
## Allocricetulus_curtatus
## Allocricetulus_eversmanni
## Alouatta_palliata
## Alouatta_pigra
## Alouatta_sara
## Alouatta_seniculus
## Alticola_albicaudus
## Alticola_argentatus
## Alticola_barakshin
## Alticola_macrotis
## Alticola_semicanus
## Alticola_strelzowi
## Alticola_tuvinicus
## Amblysomus_corriae
## Amblysomus_hottentotus
## Amblysomus_marleyi
## Amblysomus_septentrionalis
## Ametrida_centurio
## Ammospermophilus_harrisii
## Ammospermophilus_insularis
## Ammospermophilus_interpres
## Ammospermophilus_leucurus
## Ammotragus_lervia
## Amphinectomys_savamis
## Anathana_ellioti
## Andalgalomys_olrogi
## Andalgalomys_pearsoni
## Andalgalomys_roigi
## Andinomys_edax
## Anisomys_imitator
## Anomalurus_beecrofti
## Anoura_caudifer
## Anoura_cultrata
## Anoura_geoffroyi
## Anoura_latidens
## Anourosorex_yamashinai
## Antechinomys_laniger
## Antechinus_agilis
## Antechinus_bellus
## Antechinus_flavipes
## Antechinus_godmani
## Antechinus_leo
## Antechinus_minimus
## Antechinus_stuartii
## Antechinus_swainsonii
## Anthops_ornatus
## Antidorcas_marsupialis
## Antilocapra_americana
## Antilope_cervicapra
## Antrozous_dubiaquercus
## Antrozous_pallidus
## Aonyx_capensis
## Aonyx_cinerea
## Aotus_azarai
## Aotus_brumbacki
## Aotus_griseimembra
## Aotus_lemurinus
## Aotus_nancymaae
## Aotus_nigriceps
## Aotus_trivirgatus
## Aotus_vociferans
## Aplodontia_rufa
## Apodemus_alpicola
## Apodemus_epimelas
## Apodemus_gurkha
## Apodemus_hermonensis
## Apodemus_hyrcanicus
## Apodemus_iconicus
## Apodemus_ilex
## Apodemus_mystacinus
## Apodemus_pallipes
## Apodemus_ponticus
## Apodemus_semotus
## Apodemus_uralensis
## Apodemus_wardi
## Apodemus_witherbyi
## Apomys_abrae
## Apomys_camiguinensis
## Apomys_datae
## Apomys_gracilirostris
## Apomys_hylocoetes
## Apomys_insignis
## Apomys_iridensis
## Apomys_lubangensis
## Apomys_microdon
## Apomys_musculus
## Apomys_sacobianus
## Apomys_sp._SJS-2010a
## Apomys_sp._SJS-2010b
## Apomys_sp._SJS-2010c
## Apomys_sp._SJS-2010d
## Apomys_sp._SJS-2010e
## Apomys_sp._SJS-2010f
## Apomys_sp._SJS-2010g
## Aproteles_bulmerae
## Arborimus_albipes
## Arborimus_longicaudus
## Arborimus_pomo
## Archboldomys_luzonensis
## Archboldomys_sp._SAJ-2012
## Arctictis_binturong
## Arctocebus_aureus
## Arctocebus_calabarensis
## Arctocephalus_australis
## Arctocephalus_galapagoensis
## Arctocephalus_gazella
## Arctocephalus_philippii
## Arctocephalus_pusillus
## Arctocephalus_townsendi
## Arctocephalus_tropicalis
## Arctodus_simus
## Arctogalidia_trivirgata
## Arctonyx_collaris
## Ardops_nichollsi
## Arielulus_aureocollaris
## Arielulus_circumdatus
## Arielulus_cuprosus
## Ariteus_flavescens
## Artibeus_amplus
## Artibeus_anderseni
## Artibeus_aztecus
## Artibeus_bogotensis
## Artibeus_concolor
## Artibeus_fraterculus
## Artibeus_hartii
## Artibeus_hirsutus
## Artibeus_inopinatus
## Artibeus_intermedius
## Artibeus_schwartzi
## Artibeus_toltecus
## Arvicanthis_abyssinicus
## Arvicanthis_ansorgei
## Arvicanthis_nairobae
## Arvicanthis_neumanni
## Arvicanthis_niloticus
## Arvicanthis_somalicus
## Arvicola_amphibius
## Arvicola_sapidus
## Arvicola_scherman
## Asellia_tridens
## Aselliscus_stoliczkanus
## Aselliscus_tricuspidatus
## Atelerix_albiventris
## Atelerix_algirus
## Ateles_belzebuth
## Ateles_fusciceps
## Ateles_geoffroyi
## Ateles_hybridus
## Ateles_paniscus
## Atelocynus_microtis
## Atherurus_africanus
## Atherurus_macrourus
## Atilax_paludinosus
## Atlantoxerus_getulus
## Auliscomys_boliviensis
## Auliscomys_micropus
## Auliscomys_pictus
## Auliscomys_sublimis
## Avahi_cleesei
## Avahi_occidentalis
## Avahi_peyrierasi
## Avahi_unicolor
## Axis_axis
## Axis_kuhlii
## Axis_porcinus
## Babyrousa_babyrussa
## Baiomys_musculus
## Baiomys_taylori
## Balaenoptera_acutorostrata
## Balaenoptera_bonaerensis
## Balaenoptera_borealis
## Balaenoptera_brydei
## Balaenoptera_edeni
## Balaenoptera_musculus
## Balaenoptera_omurai
## Balantiopteryx_infusca
## Balantiopteryx_io
## Balantiopteryx_plicata
## Balionycteris_maculata
## Bandicota_bengalensis
## Bandicota_savilei
## Barbastella_barbastellus
## Barbastella_beijingensis
## Barbastella_darjelingensis
## Barbastella_leucomelas
## Bassaricyon_alleni
## Bassaricyon_beddardi
## Bassaricyon_gabbii
## Bassariscus_astutus
## Bassariscus_sumichrasti
## Bathyergus_janetta
## Bathyergus_suillus
## Batomys_granti
## Batomys_salomonseni
## Bdeogale_crassicauda
## Bdeogale_nigripes
## Beamys_hindei
## Beatragus_hunteri
## Belomys_pearsonii
## Berardius_arnuxii
## Berardius_bairdii
## Berylmys_berdmorei
## Berylmys_bowersi
## Bettongia_gaimardi
## Bettongia_lesueur
## Bettongia_penicillata
## Bettongia_tropica
## Bibimys_chacoensis
## Bibimys_labiosus
## Blanfordimys_afghanus
## Blanfordimys_bucharensis
## Blarina_carolinensis
## Blarina_hylophaga
## Blarinomys_breviceps
## Blastocerus_dichotomus
## Bolomys_amoenus
## Bolomys_lasiurus
## Bolomys_temchuki
## Bolomys_urichi
## Boneia_bidens
## Bos_frontalis
## Bos_gaurus
## Bos_grunniens
## Bos_indicus
## Bos_javanicus
## Bos_mutus
## Bos_primigenius
## Bos_sauveli
## Bos_taurus
## Boselaphus_tragocamelus
## Brachiones_przewalskii
## Brachylagus_idahoensis
## Brachyphylla_nana
## Brachyphylla_pumila
## Brachytarsomys_albicauda
## Brachyteles_arachnoides
## Brachyteles_hypoxanthus
## Brachyuromys_betsileoensis
## Brachyuromys_ramirohitra
## Bradypus_pygmaeus
## Bradypus_torquatus
## Bradypus_tridactylus
## Bradypus_variegatus
## Brucepattersonius_igniventris
## Brucepattersonius_soricinus
## Bubalus_bubalis
## Bubalus_carabanensis
## Bubalus_depressicornis
## Bubalus_mindorensis
## Bubalus_quarlesi
## Budorcas_taxicolor
## Bullimus_bagobus
## Bullimus_gamay
## Bullimus_luzonicus
## Bunolagus_monticularis
## Bunomys_andrewsi
## Bunomys_chrysocomus
## Burramys_parvus
## Cabassous_centralis
## Cabassous_chacoensis
## Cabassous_tatouay
## Cabassous_unicinctus
## Cacajao_ayresi
## Cacajao_hosomi
## Cacajao_melanocephalus
## Caenolestes_fuliginosus
## Calcochloris_obtusirostris
## Callicebus_coimbrai
## Callicebus_dubius
## Callicebus_nigrifrons
## Callicebus_personatus
## Callimico_goeldii
## Callithrix_argentata
## Callithrix_aurita
## Callithrix_emiliae
## Callithrix_geoffroyi
## Callithrix_humeralifera
## Callithrix_humilis
## Callithrix_jacchus
## Callithrix_kuhlii
## Callithrix_mauesi
## Callithrix_melanura
## Callithrix_penicillata
## Callithrix_pygmaea
## Callithrix_rondoni
## Callithrix_saterei
## Callorhinus_ursinus
## Callosciurus_caniceps
## Callosciurus_finlaysonii
## Callosciurus_inornatus
## Callosciurus_nigrovittatus
## Callosciurus_notatus
## Callosciurus_orestes
## Callosciurus_prevostii
## Callospermophilus_lateralis
## Callospermophilus_madrensis
## Callospermophilus_saturatus
## Calomys_callidus
## Calomys_fecundus
## Calomys_hummelincki
## Calomys_laucha
## Calomys_lepidus
## Calomys_musculinus
## Calomys_sorellus
## Calomys_sp.
## Calomys_sp._CEG40
## Calomys_tener
## Calomys_tocantinsi
## Calomys_venustus
## Calomys
## Warning in treedata(tree, data[, c("Species", "AverageGrantham", "KnKs", : The following tips were not found in 'phy' and were dropped from 'data':
## Crocidura_monticola
## Eidolon_dupreanum
## Epomops_buettikoferi
## Sorex_cylindricauda
data_tree = as.data.frame(treedata(tree_w, data[, c('Species', 'AverageGrantham', 'KnKs','DnDs' ,'GenerationLength_d')],
sort=T, warnings=T)$data)
## Warning in treedata(tree_w, data[, c("Species", "AverageGrantham", "KnKs", : The following tips were not found in 'phy' and were dropped from 'data':
## Crocidura_monticola
## Eidolon_dupreanum
## Epomops_buettikoferi
## Sorex_cylindricauda
setdiff(data$Species, data_tree$Species)# 4 вида отвалились при присоединении дерева
## [1] "Crocidura_monticola" "Eidolon_dupreanum" "Epomops_buettikoferi"
## [4] "Sorex_cylindricauda"
data_tree$AverageGrantham = as.numeric(as.character(data_tree$AverageGrantham))
data_tree$GenerationLength_d = as.numeric(as.character(data_tree$GenerationLength_d))
data_tree$KnKs = as.numeric(as.character(data_tree$KnKs))
#GenLen & Grantham
cor.test(pic(data_tree$AverageGrantham, tree_w), pic(log2(data_tree$GenerationLength_d), tree_w), method = 'spearman')
## Warning in cor.test.default(pic(data_tree$AverageGrantham, tree_w),
## pic(log2(data_tree$GenerationLength_d), : Cannot compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: pic(data_tree$AverageGrantham, tree_w) and pic(log2(data_tree$GenerationLength_d), tree_w)
## S = 6609013, p-value = 0.2717
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.05907785
cor.test(data_tree$AverageGrantham,data_tree$GenerationLength_d, method = 'spearman')
## Warning in cor.test.default(data_tree$AverageGrantham,
## data_tree$GenerationLength_d, : Cannot compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: data_tree$AverageGrantham and data_tree$GenerationLength_d
## S = 4785484, p-value = 5.302e-10
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.3245325
#GenLen & DnDs
cor.test(pic(data_tree$DnDs, tree_w), pic(log2(data_tree$GenerationLength_d), tree_w), method = 'spearman')
## Warning in cor.test.default(pic(data_tree$DnDs, tree_w),
## pic(log2(data_tree$GenerationLength_d), : Cannot compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: pic(data_tree$DnDs, tree_w) and pic(log2(data_tree$GenerationLength_d), tree_w)
## S = 6113720, p-value = 0.01556
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.1295925
cor.test(as.numeric(data_tree$DnDs),data_tree$GenerationLength_d, method = 'spearman')
## Warning in cor.test.default(as.numeric(data_tree$DnDs),
## data_tree$GenerationLength_d, : Cannot compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: as.numeric(data_tree$DnDs) and data_tree$GenerationLength_d
## S = 3984711, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.4375611
#GenLen & KnKs
cor.test(pic(data_tree$KnKs, tree_w), pic(log2(data_tree$GenerationLength_d), tree_w), method = 'spearman')
## Warning in cor.test.default(pic(data_tree$KnKs, tree_w),
## pic(log2(data_tree$GenerationLength_d), : Cannot compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: pic(data_tree$KnKs, tree_w) and pic(log2(data_tree$GenerationLength_d), tree_w)
## S = 6580489, p-value = 0.2401
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.06313875
cor.test(data_tree$KnKs,data_tree$GenerationLength_d, method = 'spearman')
## Warning in cor.test.default(data_tree$KnKs, data_tree$GenerationLength_d, :
## Cannot compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: data_tree$KnKs and data_tree$GenerationLength_d
## S = 4441184, p-value = 5.695e-13
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.3731303
######### МАНН-УИТНИ
Data = data
quantile(Data$GenerationLength_d, seq(from=0, to=1, by=0.25),na.rm = TRUE)
## 0% 25% 50% 75% 100%
## 202.795 628.295 1585.196 2726.536 18980.000
boxplot(Data[Data$GenerationLength_d < 637.3329,]$AverageGrantham, +
Data[Data$GenerationLength_d < 1825.0000 & Data$GenerationLength_d > 637.3329,]$AverageGrantham, +
Data[Data$GenerationLength_d < 2869.6933 & Data$GenerationLength_d > 1825.0000,]$AverageGrantham, +
Data[Data$GenerationLength_d > 2869.6933,]$AverageGrantham, names = c('Q1','Q2', 'Q3','Q4'), outline = FALSE, xlab = "GenLenght", ylab = "AverageGrantham",notch = TRUE)

boxplot(Data[Data$GenerationLength_d < 637.3329,]$KnKs, +
Data[Data$GenerationLength_d < 1825.0000 & Data$GenerationLength_d > 637.3329,]$KnKs, +
Data[Data$GenerationLength_d < 2869.6933 & Data$GenerationLength_d > 1825.0000,]$KnKs, +
Data[Data$GenerationLength_d > 2869.6933,]$KnKs, names = c('Q1','Q2', 'Q3','Q4'), outline = FALSE, xlab = "GenLenght", ylab = "KnKs",notch = TRUE)

##### boxplots by quartiles
boxplot(data[data$GenerationLength_d<=quantile(data$GenerationLength_d,0.25),]$AverageGrantham,
data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.25) & data$GenerationLength_d<=quantile(data$GenerationLength_d,0.5),]$AverageGrantham,
data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.5) & data$GenerationLength_d<=quantile(data$GenerationLength_d,0.75),]$AverageGrantham,
data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.75),]$AverageGrantham,
names=c('1','2','3','4'), outline = FALSE, notch = TRUE)

boxplot(data[data$GenerationLength_d<=quantile(data$GenerationLength_d,0.25),]$DnDs,
data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.25) & data$GenerationLength_d<=quantile(data$GenerationLength_d,0.5),]$DnDs,
data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.5) & data$GenerationLength_d<=quantile(data$GenerationLength_d,0.75),]$DnDs,
data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.75),]$DnDs,
names=c('1','2','3','4'), outline = FALSE, notch = TRUE)

wilcox.test(data$AverageGrantham,data$GenerationLength_d)
##
## Wilcoxon rank sum test with continuity correction
##
## data: data$AverageGrantham and data$GenerationLength_d
## W = 0, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(data$DnDs,data$GenerationLength_d)
##
## Wilcoxon rank sum test with continuity correction
##
## data: data$DnDs and data$GenerationLength_d
## W = 0, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(data$KnKs,data$GenerationLength_d)
##
## Wilcoxon rank sum test with continuity correction
##
## data: data$KnKs and data$GenerationLength_d
## W = 0, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(
data[data$GenerationLength_d<=quantile(data$GenerationLength_d,0.25),]$AverageGrantham,
data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.25) & data$GenerationLength_d<=quantile(data$GenerationLength_d,0.5),]$AverageGrantham)
##
## Wilcoxon rank sum test with continuity correction
##
## data: data[data$GenerationLength_d <= quantile(data$GenerationLength_d, 0.25), ]$AverageGrantham and data[data$GenerationLength_d > quantile(data$GenerationLength_d, 0.25) & data$GenerationLength_d <= quantile(data$GenerationLength_d, 0.5), ]$AverageGrantham
## W = 3011.5, p-value = 0.007993
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(
data[data$GenerationLength_d<=quantile(data$GenerationLength_d,0.25),]$AverageGrantham,
data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.5) & data$GenerationLength_d<=quantile(data$GenerationLength_d,0.75),]$AverageGrantham)
##
## Wilcoxon rank sum test with continuity correction
##
## data: data[data$GenerationLength_d <= quantile(data$GenerationLength_d, 0.25), ]$AverageGrantham and data[data$GenerationLength_d > quantile(data$GenerationLength_d, 0.5) & data$GenerationLength_d <= quantile(data$GenerationLength_d, 0.75), ]$AverageGrantham
## W = 2848, p-value = 0.001206
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(
data[data$GenerationLength_d<=quantile(data$GenerationLength_d,0.25),]$AverageGrantham,
data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.75),]$AverageGrantham)
##
## Wilcoxon rank sum test with continuity correction
##
## data: data[data$GenerationLength_d <= quantile(data$GenerationLength_d, 0.25), ]$AverageGrantham and data[data$GenerationLength_d > quantile(data$GenerationLength_d, 0.75), ]$AverageGrantham
## W = 2153, p-value = 1.449e-07
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(
data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.25) & data$GenerationLength_d<=quantile(data$GenerationLength_d,0.5),]$AverageGrantham,
data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.5) & data$GenerationLength_d<=quantile(data$GenerationLength_d,0.75),]$AverageGrantham)
##
## Wilcoxon rank sum test with continuity correction
##
## data: data[data$GenerationLength_d > quantile(data$GenerationLength_d, 0.25) & data$GenerationLength_d <= quantile(data$GenerationLength_d, 0.5), ]$AverageGrantham and data[data$GenerationLength_d > quantile(data$GenerationLength_d, 0.5) & data$GenerationLength_d <= quantile(data$GenerationLength_d, 0.75), ]$AverageGrantham
## W = 3478.5, p-value = 0.2975
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(
data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.25) & data$GenerationLength_d<=quantile(data$GenerationLength_d,0.5),]$AverageGrantham,
data[data$GenerationLength_d>quantile(data$GenerationLength_d,0.75),]$AverageGrantham)
##
## Wilcoxon rank sum test with continuity correction
##
## data: data[data$GenerationLength_d > quantile(data$GenerationLength_d, 0.25) & data$GenerationLength_d <= quantile(data$GenerationLength_d, 0.5), ]$AverageGrantham and data[data$GenerationLength_d > quantile(data$GenerationLength_d, 0.75), ]$AverageGrantham
## W = 2621, p-value = 0.0003169
## alternative hypothesis: true location shift is not equal to 0
hist(data$GenerationLength_d,breaks=10)

######### ЛЯМБДА
GenLen <- data_tree$GenerationLength_d
names(GenLen) <- rownames(data_tree)
Grantham <- data_tree$AverageGrantham
names(Grantham) <- rownames(data_tree)
DnDs <- as.matrix(as.numeric(data_tree$DnDs))
rownames(DnDs) <- rownames(data_tree)
phylosig(tree_w, GenLen, method = "lambda", test = TRUE)
##
## Phylogenetic signal lambda : 0.993386
## logL(lambda) : -2845.39
## LR(lambda=0) : 630.524
## P-value (based on LR test) : 3.84386e-139
phylosig(tree_w, Grantham, method = "lambda", test = TRUE)
##
## Phylogenetic signal lambda : 0.196536
## logL(lambda) : -1468.53
## LR(lambda=0) : 17.9445
## P-value (based on LR test) : 2.27441e-05
phylosig(tree_w, DnDs, method = "lambda", test = TRUE)
##
## Phylogenetic signal lambda : 0.465285
## logL(lambda) : 554.683
## LR(lambda=0) : 79.2683
## P-value (based on LR test) : 5.42225e-19
################### PGLS
MutComp = comparative.data(tree_w, data, Species, vcv=TRUE)
summary(pgls(scale(AverageGrantham) ~ log2(GenerationLength_d), MutComp, lambda="ML"))
##
## Call:
## pgls(formula = scale(AverageGrantham) ~ log2(GenerationLength_d),
## data = MutComp, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.248868 -0.046084 -0.001791 0.044802 0.212354
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.000
## lower bound : 0.000, p = 1
## upper bound : 1.000, p = < 2.22e-16
## 95.0% CI : (NA, 0.293)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.089236 0.425581 -4.9091 1.411e-06 ***
## log2(GenerationLength_d) 0.199745 0.040386 4.9459 1.183e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07277 on 347 degrees of freedom
## Multiple R-squared: 0.06585, Adjusted R-squared: 0.06316
## F-statistic: 24.46 on 1 and 347 DF, p-value: 1.183e-06
summary(pgls(scale(KnKs) ~ log2(GenerationLength_d), MutComp, lambda="ML"))
##
## Call:
## pgls(formula = scale(KnKs) ~ log2(GenerationLength_d), data = MutComp,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.32732 -0.05094 -0.00464 0.04000 0.31108
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.334
## lower bound : 0.000, p = 0.00070586
## upper bound : 1.000, p = < 2.22e-16
## 95.0% CI : (0.086, 0.610)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.708434 0.870111 -3.1127 0.0020075 **
## log2(GenerationLength_d) 0.264797 0.075496 3.5074 0.0005119 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07904 on 347 degrees of freedom
## Multiple R-squared: 0.03424, Adjusted R-squared: 0.03146
## F-statistic: 12.3 on 1 and 347 DF, p-value: 0.0005119
summary(pgls(scale(DnDs) ~ log2(GenerationLength_d), MutComp, lambda="ML"))
##
## Call:
## pgls(formula = scale(DnDs) ~ log2(GenerationLength_d), data = MutComp,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.29521 -0.03523 0.00188 0.03274 0.31114
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.331
## lower bound : 0.000, p = 3.17e-12
## upper bound : 1.000, p = < 2.22e-16
## 95.0% CI : (0.171, 0.534)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.109900 0.823284 -3.7774 0.0001864 ***
## log2(GenerationLength_d) 0.282557 0.071471 3.9535 9.337e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07493 on 347 degrees of freedom
## Multiple R-squared: 0.0431, Adjusted R-squared: 0.04034
## F-statistic: 15.63 on 1 and 347 DF, p-value: 9.337e-05
######################## ДЕЛЕНИЕ ПО СЕМЕЙСТВАМ
OrderFreq = data.frame(table(data$Order));
FrequentOrders = OrderFreq[OrderFreq$Freq >= 3,]$Var1; length(FrequentOrders) # 3!!!
## [1] 8
OrderFreq = OrderFreq[OrderFreq$Var1 %in% FrequentOrders,]
names(OrderFreq)=c('Order','NumberOfSpecies')
Mammalia = data[data$Order %in% FrequentOrders,]
agg = aggregate(list(Mammalia$AverageGrantham,Mammalia$GenerationLength_d,Mammalia$DnDs,Mammalia$KnKs), by = list(Mammalia$Order), FUN = median)
names(agg) = c('Order','AverageGrantham','GenerationLength_d','DnDs','KnKs')
test1 = cor.test(agg$AverageGrantham,agg$GenerationLength_d,method = 'spearman') ### 0.008503, Rho = 0.8434347 !!!#не работает с пропусками в продолжительности жизни
## Warning in cor.test.default(agg$AverageGrantham, agg$GenerationLength_d, :
## Cannot compute exact p-value with ties
plot(agg$GenerationLength_d,agg$AverageGrantham)

test2 = cor.test(agg$DnDs,agg$GenerationLength_d,method = 'spearman') ### 0.01071, Rho = 0.8571429 !!!
plot(agg$DnDs,agg$AverageGrantham)

test3 = cor.test(agg$KnKs,agg$GenerationLength_d,method = 'spearman') ### 0.004563, Rho = 0.9047619 !!!
plot(agg$KnKs,agg$AverageGrantham)

OrderGrantham_GenLen = c(test1$p.value,test1$estimate)
OrderDnDs_GenLen = c(test2$p.value,test2$estimate)
OrderKnKs_GenLen = c(test3$p.value,test3$estimate)
OrderFinal = data.frame(OrderGrantham_GenLen,OrderDnDs_GenLen,OrderKnKs_GenLen); rownames(OrderFinal) = c('p-value','rho')
FamFreq = data.frame(table(data$family));
FrequentFamilies = FamFreq[FamFreq$Freq >= 3,]$Var1; length(FrequentFamilies) # 3!!!
## [1] 32
FamFreq = FamFreq[FamFreq$Var1 %in% FrequentFamilies,]
names(FamFreq)=c('Family','NumberOfSpecies')
# Family
Mammalia1 = data[data$family %in% FrequentFamilies,]
agg1 = aggregate(list(Mammalia1$AverageGrantham,Mammalia1$GenerationLength_d,Mammalia1$DnDs,Mammalia1$KnKs), by = list(Mammalia1$family), FUN = median)
names(agg1) = c('Family','AverageGrantham','GenerationLength_d','DnDs','KnKs')
test4 = cor.test(agg1$AverageGrantham,agg1$GenerationLength_d,method = 'spearman') ### 0.00281, Rho = 0.5108606 !!!
plot(agg1$GenerationLength_d,agg1$AverageGrantham)

test5 = cor.test(agg1$DnDs,agg1$GenerationLength_d,method = 'spearman') ### 0.003045, Rho = 0.5128299 !!!
plot(agg1$DnDs,agg1$AverageGrantham)

test6 = cor.test(agg1$KnKs,agg1$GenerationLength_d,method = 'spearman') ### 0.01945, Rho = 0.4109991 !!!
## Warning in cor.test.default(agg1$KnKs, agg1$GenerationLength_d, method =
## "spearman"): Cannot compute exact p-value with ties
plot(agg1$KnKs,agg1$AverageGrantham)

FamilyGrantham_GenLen = c(test4$p.value,test4$estimate)
FamilyDnDs_GenLen = c(test5$p.value,test5$estimate)
FamilyKnKs_GenLen = c(test6$p.value,test6$estimate)
FamilyFinal = data.frame(FamilyGrantham_GenLen,FamilyDnDs_GenLen,FamilyKnKs_GenLen); rownames(FamilyFinal) = c('p-value','rho')
agg = merge(agg,OrderFreq)
agg = agg[order(agg$GenerationLength_d),]##"Eulipotyphla","Rodentia","Didelphimorphia","Lagomorpha","Chiroptera","Carnivora","Cetartiodactyla","Primates"
agg$AverageGrantham = round(agg$AverageGrantham,2)
agg$GenerationLength_d = round(agg$GenerationLength_d,0)
#################### Боксплоты по отрядам и семействам
Order=c("Eulipotyphla","Rodentia","Didelphimorphia","Lagomorpha","Chiroptera","Carnivora","Cetartiodactyla","Cetartiodactyla")#
median_DnDs = c(median(Data[Data$Order == 'Eulipotyphla',]$DnDs),median(Data[Data$Order == 'Rodentia',]$DnDs),median(Data[Data$Order == 'Didelphimorphia',]$DnDs),
median(Data[Data$Order == 'Lagomorpha',]$DnDs),median(Data[Data$Order == 'Chiroptera',]$DnDs),median(Data[Data$Order == 'Carnivora',]$DnDs),
median(Data[Data$Order == 'Cetartiodactyla',]$DnDs),median(Data[Data$Order == 'Cetartiodactyla',]$DnDs))
median_KnKs = c(median(Data[Data$Order == 'Eulipotyphla',]$KnKs),median(Data[Data$Order == 'Rodentia',]$KnKs),median(Data[Data$Order == 'Didelphimorphia',]$KnKs),
median(Data[Data$Order == 'Lagomorpha',]$KnKs),median(Data[Data$Order == 'Chiroptera',]$KnKs),median(Data[Data$Order == 'Carnivora',]$KnKs),
median(Data[Data$Order == 'Cetartiodactyla',]$KnKs),median(Data[Data$Order == 'Cetartiodactyla',]$KnKs))
median_GenerationLength_d = c(median(Data[Data$Order == 'Eulipotyphla',]$GenerationLength_d),median(Data[Data$Order == 'Rodentia',]$GenerationLength_d),median(Data[Data$Order == 'Didelphimorphia',]$GenerationLength_d),
median(Data[Data$Order == 'Lagomorpha',]$GenerationLength_d),median(Data[Data$Order == 'Chiroptera',]$GenerationLength_d),median(Data[Data$Order == 'Carnivora',]$GenerationLength_d),
median(Data[Data$Order == 'Cetartiodactyla',]$GenerationLength_d),median(Data[Data$Order == 'Cetartiodactyla',]$GenerationLength_d))
median_AverageGrantham = c(median(Data[Data$Order == 'Eulipotyphla',]$AverageGrantham),median(Data[Data$Order == 'Rodentia',]$AverageGrantham),median(Data[Data$Order == 'Didelphimorphia',]$AverageGrantham),
median(Data[Data$Order == 'Lagomorpha',]$AverageGrantham),median(Data[Data$Order == 'Chiroptera',]$AverageGrantham),median(Data[Data$Order == 'Carnivora',]$AverageGrantham),
median(Data[Data$Order == 'Cetartiodactyla',]$AverageGrantham),median(Data[Data$Order == 'Cetartiodactyla',]$AverageGrantham))
number_of_species = c(length(unique(Data[Data$Order == 'Eulipotyphla',]$Species)),length(unique(Data[Data$Order == 'Rodentia',]$Species)),length(unique(Data[Data$Order == 'Didelphimorphia',]$Species)),
length(unique(Data[Data$Order == 'Lagomorpha',]$Species)),length(unique(Data[Data$Order == 'Chiroptera',]$Species)),length(unique(Data[Data$Order == 'Carnivora',]$Species)),
length(unique(Data[Data$Order == 'Cetartiodactyla',]$Species)),length(unique(Data[Data$Order == 'Cetartiodactyla',]$Species)))
Genetics_Ecology = data.frame(Order,median_GenerationLength_d,median_AverageGrantham,median_DnDs,median_KnKs,number_of_species)
Genetics_Ecology = Genetics_Ecology[order(Genetics_Ecology$median_GenerationLength_d),]
All_statistics = c(GenerationLength_FractionOfSyn_fit$estimate,GenerationLength_FractionOfNonsyn_fit$estimate,GenerationLength_SummOfAllGrantham_fit$estimate,
GenerationLength_MedianOfAllNonsyn_fit$estimate,GenerationLength_AverageGrantham_fit$estimate,GenerationLength_KnKs_fit$estimate,GenerationLength_DnDs_fit$estimate)
All_pvalue = c(GenerationLength_FractionOfSyn_fit$p.value,GenerationLength_FractionOfNonsyn_fit$p.value,GenerationLength_SummOfAllGrantham_fit$p.value,
GenerationLength_MedianOfAllNonsyn_fit$p.value,GenerationLength_AverageGrantham_fit$p.value,GenerationLength_KnKs_fit$p.value,GenerationLength_DnDs_fit$p.value)
Variable = c('FractionOfSyn','FractionOfNonsyn','SummOfAllGrantham','MedianOfAllNonsyn','AverageGrantham','KnKs','DnDs')
Corr_with_GenLen = data.frame(Variable, All_statistics,All_pvalue)
boxplot(Data[Data$Order == "Eulipotyphla",]$KnKs,+
Data[Data$Order == "Rodentia",]$KnKs, +
Data[Data$Order == "Chiroptera",]$KnKs, +
Data[Data$Order == "Carnivora",]$KnKs, +
Data[Data$Order == "Primates",]$KnKs, +
Data[Data$Order =="Cetartiodactyla",]$KnKs, names = c('Eulipotyphla',"Rodentia",'Chiroptera','Carnivora',"Primates",'Cetartiodactyla'), outline = FALSE, xlab = "GenLenght", ylab = "KnKs",notch = TRUE)

boxplot(Data[Data$Order == "Eulipotyphla",]$AverageGrantham,+
Data[Data$Order == "Rodentia",]$AverageGrantham, +
Data[Data$Order == "Chiroptera",]$AverageGrantham, +
Data[Data$Order == "Carnivora",]$AverageGrantham, +
Data[Data$Order == "Primates",]$AverageGrantham, +
Data[Data$Order =="Cetartiodactyla",]$AverageGrantham, names = c('Eulipotyphla',"Rodentia",'Chiroptera','Carnivora',"Primates",'Cetartiodactyla'), outline = FALSE, xlab = "GenLenght", ylab = "AverageGrantham",notch = TRUE)

TempData = subset(Data, Data$Order == 'Eulipotyphla'|Data$Order =="Rodentia"|Data$Order =='Chiroptera'|Data$Order =='Carnivora'|Data$Order =="Primates"|Data$Order =='Cetartiodactyla')
ggplot(TempData, aes(x=Order, y=KnKs,fill=Order),notch=T) +
geom_violin(trim=FALSE)+
stat_summary(fun=median, geom="crossbar", size=0.2, color="red")+
labs(title="KnKs~GenLen",x="Order", y = "KnKs")+
theme_classic()

ggplot(TempData, aes(x=Order, y=AverageGrantham,fill=Order),notch=T) +
geom_violin(trim=FALSE)+
stat_summary(fun=median, geom="crossbar", size=0.2, color="red")+
labs(title="AverageGrantham~GenLen",x="Order", y = "AverageGrantham")+
theme_classic()

median(Data[Data$Order == "Primates",]$GenerationLength_d)
## [1] 3820.017
median(Data[Data$Order =="Cetartiodactyla",]$GenerationLength_d)# продолжительность жизни китопарнокопытных чуть больше, чем у приматов
## [1] 3595.33
# Carnivora Cetartiodactyla Chiroptera Eulipotyphla Primates Rodentia
nrow(Mammalia)
## [1] 348
for (i in 1:nrow(Mammalia)) { Mammalia$FamilyShort[i] = paste(unlist(strsplit(Mammalia$family[i],''))[c(1:3)],collapse = '') }
for (i in 1:nrow(Mammalia)) { Mammalia$OrderShort[i] = paste(unlist(strsplit(Mammalia$Order[i],''))[c(1:3)],collapse = '') }
library(boxplotdbl) # install.packages('boxplotdbl') #здесь надо разбиение по КК 0 и 1
X = data.frame(as.factor(Mammalia$FamilyShort),Mammalia$AverageGrantham)
Y = data.frame(as.factor(Mammalia$FamilyShort),Mammalia$GenerationLength_d)
par(mar = c(4, 4, 4, 4))
boxplotdou(Y,X,ylim = c(0,100), xlim = c(0,20000), name.on.axis = FALSE, cex = 1, pch = 0, cex.lab = 1, cex.axis = 1, col = rainbow(11)) # name.on.axis = FALSE factor.labels = FALSE, draw.legend = TRUE

X = data.frame(as.factor(Mammalia$FamilyShort),Mammalia$DnDs)
Y = data.frame(as.factor(Mammalia$FamilyShort),Mammalia$GenerationLength_d)
par(mar = c(4, 4, 4, 4))
boxplotdou(Y,X, xlim = c(0,15000), name.on.axis = FALSE, cex = 1, pch = 0, cex.lab = 1, cex.axis = 1, col = rainbow(11)) # name.on.axis = FALSE factor.labels = FALSE, draw.legend = TRUE

X = data.frame(as.factor(Mammalia$FamilyShort),Mammalia$KnKs)
Y = data.frame(as.factor(Mammalia$FamilyShort),Mammalia$GenerationLength_d)
par(mar = c(4, 4, 4, 4))
boxplotdou(Y,X,ylim = c(0,1), xlim = c(0,18500), name.on.axis = FALSE, cex = 1, pch = 0, cex.lab = 1, cex.axis = 1, col = rainbow(11)) # name.on.axis = FALSE factor.labels = FALSE, draw.legend = TRUE

X = data.frame(as.factor(Mammalia$OrderShort),Mammalia$AverageGrantham)
Y = data.frame(as.factor(Mammalia$OrderShort),Mammalia$GenerationLength_d)
par(mar = c(4, 4, 4, 4))
boxplotdou(Y,X, xlim = c(0,6500), name.on.axis = FALSE, cex = 1, pch = 0, cex.lab = 1, cex.axis = 1, col = rainbow(11)) # name.on.axis = FALSE factor.labels = FALSE, draw.legend = TRUE

X = data.frame(as.factor(Mammalia$OrderShort),Mammalia$DnDs)
Y = data.frame(as.factor(Mammalia$OrderShort),Mammalia$GenerationLength_d)
par(mar = c(4, 4, 4, 4))
boxplotdou(Y,X, xlim = c(0,6500), name.on.axis = FALSE, cex = 1, pch = 0, cex.lab = 1, cex.axis = 1, col = rainbow(11)) # name.on.axis = FALSE factor.labels = FALSE, draw.legend = TRUE

X = data.frame(as.factor(Mammalia$OrderShort),Mammalia$KnKs)
Y = data.frame(as.factor(Mammalia$OrderShort),Mammalia$GenerationLength_d)
par(mar = c(4, 4, 4, 4))
boxplotdou(Y,X, xlim = c(0,6500), name.on.axis = FALSE, cex = 1, pch = 0, cex.lab = 1, cex.axis = 1, col = rainbow(11)) # name.on.axis = FALSE factor.labels = FALSE, draw.legend = TRUE

##########################